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ABSTRACT 

Wc present numerical simulations of the adiabatic interaction of a shock with a clumpy 
region containing many individual clouds. Our work incorporates a sub-grid turbulence 
model which for the first time makes this investigation feasible. We vary the Mach 
number of the shock, the density contrast of the clouds, and the ratio of total cloud 
mass to inter-cloud mass within the clumpy region. Cloud material becomes incorpo- 
rated into the flow. This "mass-loading" reduces the Mach number of the shock, and 
leads to the formation of a dense shell. In cases in which the mass-loading is sufficient 
the flow slows enough that the shock degenerates into a wave. The interaction evolves 
through up to four stages: initially the shock decelerates; then its speed is nearly con- 
stant; next the shock accelerates as it leaves the clumpy region; finally it moves at a 
constant speed close to its initial speed. Turbulence is generated in the post-shock flow 
as the shock sweeps through the clumpy region. Clouds exposed to turbulence can be 
destroyed more rapidly than a similar cloud in an "isolated" environment. The lifetime 
of a downstream cloud decreases with increasing cloud-to-intercloud mass ratio. We 
briefly discuss the significance of these results for starburst superwinds and galaxy 
evolution. 

Key words: hydrodynamics - shock waves - turbulence - ISM: clouds - ISM: kine- 
matics and dynamics - galaxies: starburst. 



1 INTRODUCTION 

The energy input from high-mass stars in regions of vig- 
orous star formation creates bubbles of hot plasma which 
may overlap to create superbubbles. Such structures are 
large enough that they can burst out of their host galaxies, 
and vent mass and energy into the intergalactic medium. 
However, the properties of such flows may be dominated by 
their interaction with components with small volume filling 
fractions. A complete understanding of the gas dynamics 
of starburst regions, therefore, requires knowledge of how 
hot plasma interacts with the cold, dense molecular mate- 
rial present in the interstellar medium and in the superwind. 

Material in cold, dense clouds may be incorpo- 
rated into the hot phase via a variety of processes, in- 
cluding hydrodynamic ablation and thermal or photoion- 
ization induced evaporation. This "mass-loading" is a key 
ingre dient in models of galaxy formation and evolution 
fe.g. JSales et alj|20ich . but currently is not calculated self- 
consistently in these models. The mass-loss rates of clouds 
are also an important input parameter in models of the 
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ISM (|McKee fc Ostrikerl Il977h . Alternatively, clouds may 
be compressed by the high pressure of the hot phase and 
collapse to form stars. 

The response of a cloud to a given set of conditions re- 
mains to be fully elucidated. Although a large l iterature on 



interactions involving single clouds exists (e.g . Klein et al.l 
1 1994 lOrlando et al. 1 120081 : iPittard et all 120091 l2O10h . there 
are only a handful of stu dies on the interaction of a flow over- 
running multiple clouds ( Jun ct al. 1996; St cffen et alj|l997l ; 
iPoludnenko et alj|2002t IPittard et al.ll2005h . Thus, there re- 
mains a clear need for further studies to examine the col- 
lective effects of clouds on a flow, and to determine whether 
there are differences between the behaviours of upstream 
and downstream clouds. 

This is a first step in a long term study of the global 
effects of clouds on an overruning flow. To simplify the prob- 
lem and reduce the number of free parameters we have as- 
sumed adiabatic behaviour, with a ratio of specific heats 
7 = 5/3. It alows us to follow the generic behaviour of a flow 
interacting with clouds. In future we will perform calcula- 
tions which include heating and radiative losses. We inves- 
tigate how cloud destruction affects the density, speed and 
Mach number of the flow. We also determine whether the 
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position of a cloud within the clumpy region significantly 
affects its evolution. In Section [2] we describe the numeri- 
cal code and the initial conditions and a convergence study 
which informs the subsequent work in this paper. Our results 
are presented in Section [3l and conclusions and motivation 
for further work are addressed in Section [5] 



2 METHOD 

2.1 Numerical Setup 

The computations were performed with the mg hydrody- 
namic code. This implements adaptive mesh refinement 
(AMR), a Godunov solver and piecewise linear cell interpo- 
lation. The code solves numerically the Euler equations of in- 
viscid flow. These equations are supplemented by a sub-grid 
turbulent viscosity model, which introduces two additional 
variables; k, the turbulent kinetic energy, and e, its dissipa- 
tion rate. The fluxes of k and e across the cell boundaries are 
also returned by the Riemann solver. The values of k and e 
in each cell are updated accordingly, including changes due 
to source terms. The fu ll set of equations solved can be found 
Pittard et alj (I2009T I. The energy source term in eq. 7 of 



iPittard et alj (|2009h . Sr. is — Sk, defined in eq. 13. T is tern 
perature and r is the turbulent stress tensor. More details 
abou t the i mplementation of the k-e model can be found in 
iFalld (| 19941 ). 

The adopted grids are cartesian, and either 2-D, for 
clouds that are initially infinite cylinders, or 3-D, for clouds 
that are initially spheres. We first performed a resolution 
study for the interaction of a flow with a single cloud. Con- 
stant inflow was imposed from the negative ^-direction to- 
gether with zero-gradient conditions on the other bound- 
aries. The simulations were terminated before any of the 
cloud material reached the edges of the grid. In the multi- 
cloud simulations we imposed periodic boundary conditions 
in the {/-direction, with the same inflow and free outflow in 
the ^-direction as in the single cloud case. 

Two grids (G° and G 1 ) cover the entire domain. 
Finer grids were added where they were needed and removed 
where they were not, based on the refinement criteria con- 
trolled by differences in the solutions on the coarser grids. 
Each refinement level increased the resolution in all direc- 
tions by a factor of 2. The time-step on grid G" was Ato/2n 
where Ato was the time-step on G°. The presence of the 
k-e subgrid model adds a viscous stability condition. The 
time step is then the smaller of the viscous and hyperbolic 
time steps. An explicit method is used to advance the solu- 
tion forwards in time. In single cloud resolution tests, 4 to 8 
levels of refinement were used with the coarsest G° grid cov- 
ering the cloud with 0.5 — 2 cells per cloud radius, in order 
to achieve an effective resolution of up to 256 cells per cloud 
radius. The effective resolution is given by the resolution in 
the finest grid and is quoted as R cr , where cr is the num- 
ber of cells per cloud radius in the finest grid. For example, 
-R256 indicates an effective resolution of 256 cells per cloud 
radius. Four levels of refinement were used in the multiple 
cloud simulations to achieve a resolution of R$. All length 
scales are measured in units of a cloud radius. 

The clouds and ambient medium are initially at uni- 
form pressure. Cloud edges have a density profile from 



IPittard et ail (|2009h . with pi = 10. In the single cloud sim- 
ulations the cloud was centered on the grid origin, and a 
planar shock front was imposed at x = —3. Time zero 
is defined to be the time when the shock first encoun- 
ters the cloud. Each simulation is described by the Mach 
number of the shock, M, and the cloud density contrast, 
X- The time is measured in units of the cloud crushing 
timescale, t cc = x 1 r c i/vh, where r c i is the initial cloud 
radius and v b is t he shock velocity in the ambient medium 
(|Klein et alj 1 1994 ). Kelvin-Helmholtz (KH) and Rayleigh- 
Taylor (RT) instabilities destroy the cloud in ~ 10t cc for 
stro ng-shock interaction s, with weaker shocks taking longer 
(see IPittard et ail |2010| ) . 

In the multiple cloud simulations the clouds are 
initially randomly distributed within a rectangular region 
which we refer to as the "clumpy region" . Each cloud has an 
identical radius and density contrast, and because of the 2D 
geometry each cloud represents an infinite cylinder. The size 
of this region is 400 r-( cl ) x 400 r( c i) , unless otherwise noted. 
A key parameter in the multi-cloud simulations is the ratio 
of total cloud mass to inter-cloud mass in the clumpy region, 
MR. In this work, mass ratios of MR — 0.25 — 4 are inves- 
tigated. Time zero occurs when the shock is just outside the 
clumpy region, but is often shifted in the subsequent analy- 
sis to coincide with the shock entering the region of interest. 
The clumpy region is divided into 4 equally sized blocks 
(shown in Fig. [7} with block 1 being furthest upstream and 
the first to be hit by the shock. Various properties of the 
total cloud material within a given block and select individ- 
ual clouds within each block were monitored through the 
use of advected scalars which "colour" the flow. In pastic- 
ular, the velocity and mass of cloud material in individual 
blocks and clouds are studied. An algorithm which searches 
for the cells which are furthest downstream and contain a 
pressure higher than the ambient pressure is used to track 
the position of the shock as a function of time. 



2.2 Convergence Studies 

Previous r esolution tests of num erical shock-cloud interac- 
tions (e.g.. iNakamura et al.|[2006l ')^ have revealed that adia- 
batic, hydrodynamic simulations need ~ 100 cells per cloud 
radius (Tiioo) for a converged result. Simulations including 
more complex physics, especially ones with str ong cooling, 
have been found to need higher res olutions (|Yirak et alj 
l2010h . In contrast. Fittard et al] (|2009l ') found that adiabatic 
simulations with a k-e subgrid turbulence model show better 
convergence, and (most importantly for this work) converge 
at lower resolutions. 

Currently, multiple cloud simulations cannot be per- 
formed at a resolution as high as i?ioo and 3D simulations 
are also much more expensive computationally. The first 
(and to our knowledge, only) resolution tests of a multi- 
ple cl oud interaction were performed by iPoludnenko et aTl 
(|2004l ). who showed how resolutions of -R16, R32 and Rg4 
affect the shock position at a given time. They concluded 
that the differences between their results were smaller than 
the sensitivity of their proposed experimental design so no 
further analysis was done. All the more extensive resolution 
tests in the current astrophysical literature concern 2D sim- 
ulations of an individual cloud. We have therefore performed 
various resolution tests including 2D single cloud simulations 
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Figure 1. Convergence tests for 2D simulations of a Mach 5 shock hitting a cylindrical cloud with density contrast \ = 10 (left column), 
X = 10 2 (center column), and \ = 10 3 (right column). The time evolution of the core mass (top row) and the mean cloud velocity 
(bottom row) are shown. 



with and without a k-e model for different values of x an d 
M. A 3D single cloud simulation with M — 5 and \ — 10 
was investigated with different resolutions up to R$4 ■ Multi- 
ple (13 and 48) cloud simulations with M = 3, % = 10 2 an d 
MR = 1 and 4 were also performed at resolutions R4, Rs, 
Rie and -R32. 

We focus on two measures which are affected by the 
mixing of cloud material into the flow. These are the mean 
cloud velocity (< v x >) and the core mass of the cloud 
(Tn CO re). The latter is defined as the sum total of cloud mass 
in grid cells where more than half the cell's mass is cloud 
material. 
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2.2.1 Single cloud resolution tests 

Fig. [T] shows the evolution of m CO rc and < v x > for 2D 
simulations of cylindrical clouds hit by an M — 5 shock for 
a range of values of x- 

In all cases the interaction leads to most of the core 
mass being mixed with ambient material by t « 10t cc . The 
simulations at lower density contrasts are most resolution 
dependent. When \ = 10 the cloud mixes faster at lower 
resolutions. The opposite is true when x = 10 2 - For x i$ 10 2 , 
simulations with at least 32 cells per cloud radius are much 
better converged than simulations at lower resolution. For 
the highest density contrast studied (x = 10 3 ) convergence 
is obtained at lower resolution. 

We now examine the resolution dependence of a 
spherical cloud in a 3D simulation. Fig. [2] shows the evo- 
lution of the core mass and mean cloud velocity as a func- 
tion of resolution for a cloud with x = 10 2 hit by a Mach 5 
shock. While the velocity diverges a bit after 5t cc , the core 
mass evolution at R\§ exhibits the same features as at R32 
and i?64, whereas in 2D a resolution of R32 is required for a 
similar level of convergence. Thus, it appears that the extra 
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Figure 2. Convergence tests for 3D simulations of a Mach 5 
shock hitting a spherical cloud with density contrast \ = 10 2 . 
The time evolution of the core mass (top) and the mean cloud 
velocity (bottom) arc shown. 



degrees of freedom for instabilities in 3D versus 2D has the 
effect that convergence is achieved at lower resolutions. 

The geometry makes a significant difference to the 
cloud evolution. Panel b of Fig. [1] shows that a cylindri- 
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Figure 3. Comparison of the core mass evolution when using the 
linear/exact Riemann solver in mg and other Riemann solvers. 
The k-e subgrid model is enabled in the top panel and disabled 
in the bottom panel. We experienced difficulties with one of the 
Tadmor calculations which stops at t = 3 t cc . 



cal cloud has lost about half of its core mass by t ~ 8t cc , 
whereas Fig. [2] shows that for a spherical cloud this occurs 
by t w 6 t cc ■ 

If the time-averaged pressure returned by the linear 
solver differs by less than 10% from the pressures in the 
left and right states, the solution from the linear solver is 
used. Otherwise an exact solver with the standard Secant 
method is used. We have compared the effects of differ- 
ent Riemann solvers on singl e cloud simula tions. In par- 
ticular we have used HLL E jEinfeldtl Il988h and Tadmor 
ijNessvahu fc Tadmorlll99oh solvers. With a first order (spa- 
tial and temporal) update the scheme is very diffusive ir- 
respective of which Riemann solver is used (although the 
mg solver is least diffusive) . Fig. [3] shows results with the 
standard second order scheme. It reveals that the different 
solvers produce similar results at high resolutions, while the 
k-e model red uces the resolution dependance, especially at 
late times (cf. iPittard et ai]|2009l '). 



2.2.2 Multiple clouds resolution test 

We have also performed a convergence study for 2D multi- 
cloud simulations. In this study a Mach 3 shock overruns 48 
identical cylindrical clouds, each with a density contrast to 
the intercloud medium of \ = 10 2 - The ratio of cloud mass 
to intercloud mass within the clumpy region is MR — 4. 
Fig.Ushows the initial setup and distribution of clouds. The 
clouds fill a region which is 100 r c \ deep and 40 r c i wide, with 
periodic boundaries in y. Resolutions of R4, R$, R\q and -R32 



Figure 4. The initial (t = 0) cloud distribution for a simulation of 
a shock overrunning multiple cylindrical clouds. The shock Mach 
number M = 3, the cloud density contrast \ = 10 2 , and the 
ratio of cloud mass to intercloud mass in the clumpy region is 
MR = 4. The arrow marks an individual cloud for which the 
properties were monitored in time - see Fig. [6] In this plot and 
in all other density snapshots the spatial scale is in the units of 
the cloud radius and the density is in the units of the ambient 
density. 



were used. The latter required a grid with an effective size 
of 12800 x 1280. 

The shock exits the clumpy region at around t — 
12.5 t cc in all of the simulations. Density plots for times 
just after this moment are shown in Fig. [5] for the Rg and 
R32 simulations. Higher compressions and greater fragmen- 
tation are seen in the higher resolution simulation. However, 
while significant differences in the behaviour of any individ- 
ual cloud in the simulation can be identified, Fig. [6] shows 
that important global parameters in the multiple cloud sim- 
ulation, such as the rate at which cloud mass is mixed into 
the flow (shown as the time evolution of the total mass of all 
of the cloud cores in Fig. ®l) , the rate at which momentum 
is transferred from the flow to the clouds (shown as the ve- 
locity of a single cloud, Fig.[6jj), and shock velocity (Fig.[6p) 
are not very sensitive to the resolution used. 

Indeed, it is clear from Fig.[S]that a resolution of R$ is 
sufficient in order to obtain an accurate representation of the 
global effect of multiple clouds on a flow. In previous multi- 
cloud si mulations resolutions of R 32 an d Riq have been 
adopt ed (|Poludnenko et ail (|2002T l and iPoludnenko et all 
l|2004h . respectively), but we show here that a somewhat 
lower resolution can be safely used, at least when a sub- 
grid turbulence model is employed. Therefore, we perform 
all other simulations in this work at a resolution of Rs- 
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Figure 5. Resolution test for a Mach 3 shock overrunning mul- 
tiple cylindrical clouds with the initial setup shown in Fig. [4] A 
logarithmic density plot is shown at t = 14.2 t cc for resolutions of 
R s (left) and R32 (right). 



Table 1. Summary of the multi-cloud simulations performed. 
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Figure 6. Resolution test for a Mach 3 shock overrunning mul- 
tiple cylindrical clouds with \ = 10 2 and MR = 4 (see Fig. [SJ. 
The time evolution of a) the total mass of all of the cloud cores, 
b) the mean shock velocity (normalized to the ambient intercloud 
sound speed) and c) the average core velocity of a single cloud, 
are shown. 



3 RESULTS 

In this section we show the results of a number of simulations 
with different values of M, x, an d MR. Table [1] summarizes 
the simulations that were performed. We adopt a naming 
convention such that m3chi2MRl refers to a simulation with 
M = 3, x = 10 2 and MR = 1. 



1 Number of clouds in the clumpy region. 

2 Simulation chi2MRl_double has the same cloud distribution as 
chi2MRl but it is repeated in the x-direction to obtain a distri- 
bution with twice the depth along the direction of shock propa- 
gation. 



3.1 Mach 3 shock interactions 

We first focus on simulations for a Mach 3 shock. Figs. [Jj 
and [8] show the time evolution of the density distribution 
in simulation m3chi2MR4, which contains 1959 clouds. The 
shock sweeps through the clumpy region, propagating fastest 
through the channels between clouds. Behind the shock, the 
overrun clouds are in various stages of destruction, with the 
clouds closest to the shock being those which are in the 
earliest stage of interaction and which, therefore, are the 
most intact. Further behind the shock the clouds gradually 
lose their identities as they are mixed into the flow. Fig. [7] 
shows that a global bowshock moves upstream into the post- 
shock flow. Although the bowshock disappears out of view 
in Fig. [8] it remains on the grid which extends to x = —450. 
Low density (and pressure) regions in the post-shock flow 
are visible behind clouds as the flow rushes pass. The global 
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shock front is momentarily deformed by each cloud that it 
encounters as it sweeps through the clumpy region. These 
local deformaties in the shock front gradually accumulate 
into a distortion of the shock front on larger length scales. 
The shock is most deformed as it reaches the end of the 
clumpy region (see the middle panel of Fig. [8]). The other 
striking feature of the interaction is the formation of a dense 
shell in the post-shock flow due to the addition of mass from 
the destruction of the clouds. This is clearly seen in the pan- 
els in Fig. [8] and is the most important large-scale structure 
formed in the global flow as a result of the interaction. 

The shock slows as it sweeps through the clumpy re- 
gion and transfers momentum to cloud material. When the 
clumpy region is highly porous, the deceleration of the shock 
front is minimal and the clouds gradually reach the veloc- 
ity of the postshock flow. This occurs if the global cloud 
cross section is small (for instance, when the mass ratio is 
low (e.g., MR = 0.25) and/or when the cloud density con- 
trast x i s high ( e -g-> X = 10 3 ) MR = 1 as in simulation 
mSchiSMRl). 

Fig. shows snapshots of the density distribution 
at the time that the shock exits the clumpy region in 
models m3chi2MRl (top), m3chi2MRl_double (middle), and 
m3chi3MR4 (bottom). All of these models have lower num- 
ber densities of clouds than model m3chi2MR4: in the first 
two models it is because the value of MR is reduced, while in 
model m3chi3MR4 it is because the density contrast of the 
clouds has increased while the value of MR is unchanged. 
As such, the clumpy region in each of these models is more 
porous than in model m3chi2MR4, and the shock is able 
to sweep through without such a significant reduction in its 
velocity. Comparison of Figs. [8] and [9] also reveals signif- 
icant differences between these models in the structures of 
the flow through the clumpy region. 

3.1.1 Velocity behaviour 

As the shock sweeps through a clumpy region, it causes 
clouds in block 1 to accelerate first; then clouds in blocks 2, 
3 and 4 accelerate. This transfer of momentum to the clouds 
inevitably causes the shock and the post-shock flow to slow. 
Figs. [TU] and [TJJ show the time evolution of the average ve- 
locity of cloud material within particular blocks of clouds 
(e.g., as delineated in the top panel of Fig. [7]) normalized to 
the intercloud sound speed. 

Fig. [10] shows that the shock front slows down to 
about 50% of its initial velocity (though it remains super- 
sonic). In this simulation we find that the first 3 blocks of 
clouds accelerate up to 1.25 ci (where ci is the sound speed 
in the intercloud ambient medium). However, by this point 
the post-shock flow immediately after the shock has deceler- 
ated to marginally below 1 ci . Hence, material stripped from 
the upstream clouds pushes against the downstream clouds, 
compressing the clumpy region. 

The shock front advances at roughly constant veloc- 
ity after having reached a minimum while in the clumpy re- 
gion. However, as the shock front leaves the clumpy region 
it reaccelerates at a roughly constant rate until it reaches 
a velocity somewhat smaller than the initial shock veloc- 
ity. After that it accelerates very slowly, and asymptotically 
reaches its initial velocity (corresponding to Mach 3). 

The final block of clouds is not held back by further 



clouds downstream and therefore does not stop accelerating 
at the same velocity as the other blocks. Instead it acceler- 
ates until it reaches the velocity of the post-shock flow after 
the shock stops reaccelerating. Blocks 3, 2 and 1 repeat this 
behaviour in that order. 

Figs. Illb ) and b) show the velocities obtained in 
simulations m3chi2MRl and m3chi2MRl_double, for which 
the mass ratio of cloud to intercloud material in the 
clumpy region is unity (the number density of clouds is 
therefore 4x lower than in model m3chi2MR4). Model 
rn3chi2MRl -double has a cloud distribution that is identi- 
cal to that of model m3chi2MRl, but the distribution is 
repeated once to create a clumpy region that is twice as 
deep. In both cases the clumpy region is more porous and 
therefore, the reductions in the shock and post-shock veloc- 
ities are not as severe as those seen in model m3chi2MR4- 
For this reason, the clouds are also accelerated initially to a 
higher velocity (~ 1.7ci) than in model m3chi2MR4- 

As previously stated, the clumpy region can also be 
made more porous by increasing \ f° r a fixed MR. Fig. Illb ) 
shows the velocity profiles obtained from model m3chi3MR4- 
As expected, the shock does not decelerate as much as in 
model m3chi2MR4- However, model m3chi3MR4 behaves 
differently in another way: because each cloud is more resis- 
tant to the flow, the clouds are not completely destroyed by 
the time the shock leaves the clumpy region (at t « 13t cc ) 
when the collective cloud material from each block is still 
accelerating. In turn, due to continuing mass-loading of the 
flow, the shock continues to decelerate after it has left the 
clumpy region. In fact, the shock front decelerates until 
t « 36 t cc , remains at constant velocity until t ~ 44t cc and 
reaccelerates thereafter. 

3.1.2 Stages of interaction 

The behaviour of the simulations noted in the previous sec- 
tion allows the identification of 4 distinct phases in the evo- 
lution of the shock front. Firstly, there is a shock deceleration 
phase as the shock enters the clumpy region. This may be 
followed by a steady phase, when the shock front moves at 
constant velocity. The third stage is a reacceleration phase 
during which the shock front accelerates (at a fairly steady 
rate) into the homogeneous ambient medium. This stage al- 
ways starts after the shock has traversed through the whole 
of the cloudy region, but not necessarily immediately after. 
The final stage begins once the shock's acceleration slows. 
At this point the shock propagates with a velocity slightly 
slower than its initial velocity, but is continuing to acceler- 
ate very slowly, due to the constant velocity inflow trying to 
return the shock to its initial velocity. 

The number and timings of the stages is model spe- 
cific. In models for which the lifetime of the clouds is less 
than the crossing time of the shock through the clumpy 
region, the steady phase ends when the shock leaves the 
clumpy region (see, e.g., simulation m3chi2MR4 in Fig. 1 10[) . 
In models for which significant mass loading of the flow con- 
tinues after the shock leaves the clumpy region the end of the 
steady phase is delayed (see, e.g., m3chi3MR4 in Fig. I lib ). 

The steady stage is best seen in models m3chi2MR4 
and m3chi2MRl_double, and to some extent it is also visible 
in model m3chi3MR4- These are the simulations with the 
most mass in the clouds. Models chi2MR4 and chi3MR4 
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Figure 7. The time evolution of the logarithmic density for model m3chi2MR4, shown at t = (top), t = 8.1 t cc (middle) and t = 16.2 t cc 
(bottom). The top panel also shows how the different cloud regions are defined. 



have the highest cloud to intercloud mass ratios of 4. In 
fact, the onset of a steady stage seems to intimately depend 
on the formation of a dense shell (see Sec. 13.1.3]) . 

In contrast, model m3chi2MRl (Fig. I lib ,) does not 
achieve a steady state. Instead the shock velocity evolves 
from the deceleration phase immediately into the reacceler- 
ation phase. One might expect that there is a better chance 
for a steady stage to occur if the clumpy region is deep. 



Fig. II lb indeed shows that a steady state occurs in model 
m3chi2MRl -double where the clumpy region is twice as deep 
as in model m3chi2MRl. The onset of the steady stage oc- 
curs when the shock is approximately halfway through the 
clumpy region, and all 4 stages can now be distinguished as 
in model m3chi2MR4. 

For all of the x = 10 2 models, the reacceleration stage 
begins the moment the shock leaves the clumpy region, but 
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Figure 8. As Fig. but at t = 36.6t C c (top), t = 61.0fcc (middle) and t = 85.5t C c (bottom). 



for the x — 10 runs the onset of this stage shows more 
complicated behaviour. Because clouds with a higher den- 
sity contrast evolve more slowly, significant evolution of the 
global flow can occur after the shock front leaves the clumpy 
region. Fig. lllb shows that this occurs in model m3chi3MR4, 
as the continued deceleration of the shock is evident. At 
some later time, presumably when the clouds are finally fully 
entrained into the flow, the shock starts to accelerate. There 
is also a hint of a steady stage in model m3chi3MR4, but 



examination of the velocity graphs reveals that it occurs at 
a much later evolutionary stage than in model m3chi2MR4 
(as seen, e.g., when compared to the velocity evolution of 
material mixed in from the various blocks) . Given that there 
are about 10 x fewer clouds in model m3chi3MR4 than in 
m3chi2MR4 this can be explained by the relative timescales 
involved in the shell formation which is discussed next. 



Simulations of shocks encountering clumpy regions 9 




-50 . - 



-100.0 b. 



-200 . 



. 



200 .0 




log(RHO) 



2.0 



. 



200 . 



400 .0 



600 .0 




-200. 



. 



200 .0 



log(RHO) 



3.( 



. i 



Figure 9. Snapshots of logarithmic density for models mSchi2MRl (top), m3chi2MRl -double (middle), and m3chi3MR4 (bottom). 
The snapshots are taken as the shock front is exiting the cloudy region which occurs at t = 46 tcc, t = 95 tec and t = 13.5 t cc respectively. 
Note the different colour scale of the bottom panel. 



3.1.3 Shell formation and evolution 

Figs. [7] and [8] show that in model m3chi2MR4, a shell forms 
in the post-shock flow at t = 16.2 t cc . It is fully developed 
by 36.6 t cc . At t = 36.6 tcc 7 distinct regions can be distin- 
guished. A uniform ambient medium (without any clouds) 
lies at the right edge of the figure, while at the far left of 
the grid (off the figure in this case) lies the original post- 
shock flow specified by inflow boundary conditions at the 



upstream edge. A further 5 regions (from right to left) ex- 
ist inbetween these other two. A region of unshocked clouds 
lies in the range 50 < x < 200. Shocked clouds embedded 
in a relatively low density postshock flow lie in the region 
30 < x < 50. Shocked clouds are being overrun by a denser 
shell in the region —20 < x < 30, though individual cores 
are still visible. This shell becomes gradually more uniform 
towards x ~ —70. The upstream edge of this dense shell 
is fairly distinct from the less dense gas that it abuts (at 
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Figure 10. The time evolution of the mean ^-velocity for cloud 
material for each different block in simulation m3chi2MR4- Also 
shown are the average velocity of the shock front and the velocity 
immediately behind the shock front. The vertical line indicates 
the time when the shock leaves the clumpy region, while the hor- 
izontal line indicates the initial speed of the post-shock flow. 
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Figure 11. As Fig. [10] but for simulations a) m3chi2MRl, 
m3chi2MRl_double and c) m3chi3MR4. 



b) 




Figure 12. The density, volume averaged over y, as a function of 
x in simulation M3chi2MR4, at different times. The shock exits 
the clumpy region at t si 60 t C c • 



x ~ —70). The latter is gas in the original supersonic post- 
shock flow which has shocked against the clumpy region. 
This final region is bounded by a bow shock which lies off 
the left edge of the figure. 

Fig. ll2l shows the evolution of the density in the shell. 
Initially the average postshock density in the clumpy re- 
gion is only about twice that of the shocked intercloud gas. 
However, as the clouds are destroyed and their mass mixes 
into the global flow, the density in this region increases. In 
model m3chi2MR4, a dense region, downstream of the re- 
gion where distinct clouds are still visible, has formed by 
about the time that the shock reaches its minimum velocity. 
This is the "shell". The compression of the shell progresses 
a little further, until a state in which the maximum den- 
sity is steady in time is reached. This is reinforced by the 
simulation (model m3chi2MRl -double) for a wide clumpy 
region, which shows that the shell only becomes wider with 
time after this point. Once there are no more clouds block- 
ing the path of the shell, it starts to expand from its leading 
edge, and the maximum density within it drops. This is seen 
clearly from Fig. 1121 and is preceded by the shock reaching 
the edge of the clumpy region. Simulations with \ = 10 3 
behave slightly differently - in these the density of the shell 
continues to increase after the shock has left the clumpy 
region due to the continued ablation of material from the 
longer-lived clouds. 

3.1.4 Mach number profile 

Fig.[l3]shows the y-averaged Mach number profile of the flow 
in simulation M3chi2MR4 at t = 44.8 t cc . The undisturbed 
post-shock flow is mildly supersonic (M w 1-04), as seen 
at the far left of the panel. The shock at this time is just 
downstream of x — 100. 

In the region between 70 < x < 108 clouds are 
overrun by the flow and accelerated by the gas streaming 
past, so that the Mach number of the flow increases as x 
decreases. While the clouds remain identifiable as distinct 
entities, their interiors remain colder than the shocked in- 
tercloud flow, as initially the dense cloud is in pressure equi- 
librium with its less dense surroundings, so that its tempera- 
ture is lower than that of the intercloud medium by a factor 
of x- Together these effects ensure that the Mach number of 
the flow eventually exceeds the Mach number in the undis- 
turbed post-shock flow. The local Mach number reaches a 
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Figure 13. The local Mach number, volume averaged for all 
values of y, as a function of x in simulation M3chi2MR4, at 
t = 44.8 t C c- The shock is at x ?s 108 and the bow shock is at 
x £3 -265. 



peak of about 1.4, at a location corresponding to the down- 
stream edge of the shell (which occupies a region between 
-50 < x < 70). 

As the clouds begin to lose their individual identities 
their material is heated as it is mixed into the surrounding 
flow. This increases the sound speed of this material, and re- 
duces the overall Mach number of the flow. The Mach num- 
ber continues to decline with decreasing x until the cloud 
material is fully mixed (which occurs at the upstream edge 
of the shell in this simulation). The region of constant Mach 
number (M w 0.6) between the upstream edge of the shell 
(at x w —50) and the bow shock at the head of the clumpy 
region (at x ~ —265) is the gas in the original post-shock 
flow which has shocked against the clumpy region. 

3.2 Dependence on shock Mach number 

Fig. ll4l shows density snapshots when the shock front is exit- 
ing the clumpy region for models ml.5chi2MR4, m2chi2MR4 
and ml0chi2MR4- It is clear that a slower shock produces 
less compression, and that the shell which then forms is 
wider. In addition, weak shocks take longer to destroy 
clouds, so the shell forms further from the shock front. Per- 
haps most significantly, as the shock decelerates during its 
passage through the clumpy region it can slow so much that 
it decays into a wave which advances at the intercloud sound 
speed. This is seen in both the M = 1.5 and M = 2 simula- 
tions for which Fig. [14] reveals very weak density jumps at 
the edges of the clumpy regions. 

Fig. [TS] displays the x component of the velocity av- 
eraged over all y for the ml.5chi2MR4 and m2chi2MR4 sim- 
ulations normalized to the ambient intercloud sound speed. 
The points show the velocity of the shock calculated by 
an algorithm which searches for a pressure jump. However, 
when the shock decays into a wave the pressure jump largely 
disappears and this calculation fails. The points with error 
bars then show the velocity of the disturbance determined 
by measuring, by eye, the position of the leading edge of the 
disturbance. In the M — 2 simulation the average velocity 
of the shock/wave disturbance is close to transonic. Close 
examination reveals that the shock appears fragmented and 
local regions of sonic waves can be seen. The shock com- 
pletely decays into a sonic wave in the M — 1.5 simulation. 
These results are in harmony with theoretical predictions 
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Figure 15. The time evolution of the x component of the velocity 
averaged over all values of y, for different regions in simulations 
ml.5chi2MR4 (top) and m2chi2MR4 (bottom). Also shown are 
the average shock front velocities. The shock velocity varies when 
parts of the shock front become subsonic. Vertical lines indicate 
the time when the shock leaves the clumpy region. 



based on a uniform density region suggesting that a super- 
sonic to subsonic transition occurs at this mass ratio if the 
initial shock Mach number is M ~ 2 (see Sec. l3.3l for further 
details). 

3.3 The reduction in the shock speed 

Fig- HE] shows the minimum shock speed (normalized to the 
initial shock speed) as a function of shock Mach number 
obtained in simulations with MR = 4 and x — 10 (plotted 
as black crosses). The minimum shock speed in each case 
occurs during a "steady" stage when the shock is moving 
through the clumpy region. The error bars correspond to the 
lcr spread in measured speeds. We have already noted that in 
the M = 1.5 simulation the shock slows to the point that the 
disturbance becomes a wave. In this case the minimum speed 
is limited by the sound speed in the intercloud medium, as 
seen in Fig. 1161 The situation in the M = 2 simulation is 
not so straightforward - a global shock front is visible but 
there are local regions where a wave is seen instead. 

We can compare the measured shock speed to the 
expected speed of a shock transmitted into a r egion of en- 
hance d density. For this we make use of eq. 5.4 in lKlein et alj 
(1994), for a shock interacting with a single cloud: 

= (F st F cl ) 1/2 X - 1/2 v . (1) 

Here P st = P3/P1, where Pi is the pressure far upstream and 
P3 is the stagnation pressure (just upstream of the cloud), 
and Pci = P4/P3 is the ratio of the pressure just behind the 
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Figure 14. Snapshot of logarithmic density for models ml.5chi2MR4 (top), m2chi2MR4 (middle) and ml0chi2MR4 (bottom). The 
snapshots are taken as the shock front is exiting the cloudy region which occurs at t = 53t cc , t = 62t cc and t = 58t cc respectively. 
Relative to its initial speed, the M = 2 shock is slowed the most. 



transmitted shock and the stagnation pressure. In the limit 
of high Mach number and high density contrast, 

Fst ^ 1+ l + lttx-W (2) 

To apply this expression to our multi-clump simulations we 
need an "effective" value for \ which accounts for the in- 
homogeneity of the clumpy region. We choose this value to 
be Xeff = MR + 1, which is the average density contrast of 
the region (total mass of the region divided by the mass a 



same size region of ambient density would have). Further- 
more, in simulations where the clumpy region is replaced 
with a region of uniform density (see Sec. 13. 4[) , it is clear 
that Pi = P3 = P2, where P2 is the pressure just after the 
bow shock. We shall assume that this is also the case for our 
multiple cloud simulations, and therefore adopt _F c i = 1. 
For our clumpy simulations with MR = 4 and \ = 10 2 we 
therefore obtain XsS = 5 and F st = 1.55, so that Eq.[T]gives 
v s — 0.56vo- As can be seen, this is comparable to the values 
shown in Fig. 1161 in which the v s = 0.56wo line is plotted. 
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Figure 16. The reduction in shock speed (black crosses, as mea- 
sured during the "steady" phase as the shock propagates through 
the clumpy region) normalized to the initial shock speed. The re- 
duction is plotted as a function of the initial shock Mach number, 
M, for simulations with \ = 10 2 an d MR = 4. Also shown is the 
equivalent reduction in speed of a shock propagating through a 
uniform density enhancement (pink crosses - see Sec. 13. il l. 



3.4 Comparisons to a shock encountering a region 
of uniform density 

In the limit of an infinitely deep/thick cloud distribution we 
expect the global behaviour of the shock to approach that of 
a shock encountering a uniform medium of the same average 
density. This could be similar to the steady state reached in 
some of our simulations. Therefore, we have performed a 
number of additional calculations of a shock encountering 
a uniform region of enhanced density equal to the average 
density of a clumpy region with MR — 4. The minimum 
shock speed through this region is shown in Fig.[l6]with pink 
crosses. As can be seen, the agreement with the multi-cloud 
simulations is very good, with significant deviation only at 
M = 1.5. Clearly the minimum shock speed obtained in our 
multi-cloud simulations is close to the shock speed occuring 
in a comparable region of uniform density. 

Fig. ll7l compares the x-dependance of the density, av- 
eraged over all values of y, given by simulations of a shock 
interacting with a uniform region of enhanced density to 
analogous profiles obtained from our multiple cloud simu- 
lations. One notices the similarity between the profiles: the 
maximum density of the shell and the reduction in density 
after the shock exits the higher density /clumpy region are 
similar. Obviously, the uniform region simulations yield dis- 
tinct edges and slopes corresponding to the shocks and rar- 
efaction waves, whereas the multiple cloud simulations give 
smoother profiles. The smoothing length is of the order of 
Leo, which is the distance th at a cloud is displaced b y the 
flow prior to its destruction (jPoludnenko et al"]|2002f ). For 
the same reason the width of the densest part of the shell is 
narrower in the multiple cloud simulations. 

The picture of a shock transmitted into a uniform 
density region is different. The velocity jumps according to 
Rankine-Hugoniot conditions. As the shock front is moving 
faster than the postshock flow, the width, w of the com- 
pressed region is increasing at a rate given by: 

dw 

~d~t =Vs ~ Vps 




3 Mu 
4M„ + 4 



(3) 



where M u is the Mach number of a shock transmitted into a 
uniform medium and c u is the sound speed in that medium. 



Figure 17. Snapshots of the x-dependance of the density av- 
eraged over all values of y (thick lines) from simulations a) 
m2chi2MR4, b) m3chi2MR4 and c) ml0chi2MR4. Also shown 
are density profiles from simulations with uniform but corre- 
sponding density. In all cases the shock exits the region of en- 
hanced density (smooth or clumpy) at t m 60 t C c ■ 



Eqs.[T]and[2]can be used to determine M u , but the approx- 
imation in Eq. [2] is only valid at high Mach numbers and 
high density contrasts. However, we can obtain an estimate 
based on Eq.[T] As F 3t ^ 1 and, because the shock is steady, 
F c = 1 (|Klein et al.|[l993 ). M u > M and so ^ > ^p. 
Furthermore, two limiting cases can be found. In the limit 
of a weak shock ^ ~ c u ; this is the minimum value, and it 
would be higher in all other cases. In the limit of a strong 
shock and high x, ~ 0.44M c„. 

Fig. 1171 shows that the widths of shells in clumpy re- 
gions correspond well to those in regions of equivalent uni- 
form density. As such Eq.[3]can be used, with M u either de- 
termined from a corresponding uniform density case, or by 
replacing it with v a /c u , with c u determined for the medium 
of corresponding uniform density. The limiting cases would 
be different though. In the weak shock limit, the shock may 
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decay into a sound wave. In such a case ^ «Ci, but it could 
be significantly less depending on how this affects individual 
cloud lifetimes. In the limit of very high \ the upstream edge 
of the shell takes a long time to reach v ps and is effectively 
stationary, so = v s . 

3.5 Individual cloud evolution 

We expect individual clouds to evolve differently when over- 
run by a mass-loaded evolved shock which has been altered 
by the presence of clouds further upstream than when they 
are overrun by a steady planar shock. Section 13.1.11 shows 
that the shock Mach number and in turn the velocity of 
the post shock flow decreases as the shock sweeps up cloud 
material. The lower velocity of the shock should prolong 
the lifetime of a cloud, but this may be offset by the in- 
creased density of the flow. Furthermore, the postshock flow 
is much more turbulent. Thi s aids the developme nt of cloud- 
destroying instabilities (see iPittard et al. 21)11!):. Finally, a 
region of higher density follows some distance behind the 
shock. In the extreme case it is a dense shell, but even if 
a shell does not fully form, remnants of upstream clouds 
may well interact with clouds further downstream. There- 
fore, it is difficult to predict whether downstream clouds are 
destroyed more readily than their upstream counterparts. 
However, we can use our simulation results to investigate 
this. 

Mass loss rates from different clouds in the simula- 
tions are shown in Fig. [18] Each cloud is from a particu- 
lar block (see, e.g., Fig. [7]), and the time in each plot has 
been shifted so that t = corresponds to the time when 
the shock first encounters the cloud under consideration. 
The greatest differences between the evolution of different 
clouds in a given simulation occur in model ml.5chi2MR4, 
as shown in Fig. 118b ). Here the cloud in the block which 
is first to be hit by the shock (labelled singiel) encounters 
a flow which is relatively unmodified by the small number 
of clouds which lie further upstream. Although the interac- 
tion is relatively weak because of the modest shock Mach 
number, significant mass-loss from the cloud occurs after 
t — 5t cc , and the core of the cloud is completely destroyed 
by t m 27 fee (a similar lifetime occurs for spherical clouds 
- see lPittard et akll2010h . In contrast, significant mass-loss 
does not occur until t » 10 f cc for the cloud in the second 
block (labelled single2), and until t ~ 12i cc and 13 £ cc for 
clouds in the third and fourth blocks (labelled single3 and 
single4, respectively). The low initial mass- loss rates from 
clouds further downstream is due to the decay of the shock 
into a wave as it moves through the clumpy medium. While 
the flow past the cloud is subsonic the mass-loss rate remains 
very low. 

The same initial behaviour is also seen in simulation 
m2chi2MR4, as shown in Fig. 118b ). In contrast, at higher 
shock Mach numbers (models m3chi2MR4 and ml0chi2MR4 
in Fig. I18p the shock remains supersonic as it transits 
through the clumpy region, and differences between the 
clouds in the evolution of their mass-loss rates are not so 
readily apparent (see Fig. I18fc and d) . 

It is clear that the cloud lifetimes must be compared 
in a statistical way. Fig.[T5]shows the ratio of the lifetimes of 
individual clouds to the lifetime of the equivalent cloud hit 
by a "clean" shock (i.e. initially with no post-shock struc- 
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Figure 19. Ratio of the lifetime of a downstream cloud to the 
equivalent lifetime obtained when the cloud is hit by a "clean" 
shock with a uniform post-shock flow. 



ture due to interactions with upstream clouds), 7iiif c . The 
top panel in Fig. [T5] shows the results for M — 3, while 
the bottom panel shows the results for M = 10. In both 
cases there is a general trend for a reduction in Rm c with 
increasing MR, albeit with a reasonably large scatter for 
a particular cloud. Table [2] notes the average value of R\a e 
in a given simulation, where this trend is now even clearer. 
When MR — 4, we find that downstream clouds have lives 
which are about 40% shorter than the lifetime of an iso- 
lated. Table [2] also indicates that there is no clear trend 
of Rnia with the shock Mach number M for a given value 
of the mass-ratio, MR. Finally, a simulation performed at 
twice our standard resolution reveals that there is currently 
a slight resolution dependence on i?iif c (model m3chi2MR4 
returns R n{e = 0.62 ± 0.07 and 0.72 ± 0.02 at a resolution of 
8 and 16 cells per cloud radius, respectively). 

The reduction in the cloud lifetime due to the tur- 
bulent nature of the flow is clearly more significant than 
the weakening and broade ning of the shock which acts to 
increase the cloud lifetime (| Williams fc Dvson|[2002T ). 



4 DISCUSSION 

Our work is relevant to objects where hot, diffuse gas in- 
teracts with a colder dense phase. In many of these objects 
the cold phase, despite its low volume filling fraction, may 
dominate the dynamics of the hot phase, and thus signifi- 
cantly change the morphology and evolution of the object, 
and its emission. On the smallest scales these ob jects include 
wind- blown bubbles and supernova remnants. ICowie et al.1 
(|l98ll ) were the first to study the behaviour of supernova 
remnants expanding into a clumpy medium. They found 
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Table 2. Ratio of the lifetimes of a downstream cloud to the 
equivalent lifetime of an "isolated" cloud. The quoted ratio is 
averaged over the lifetime of the individual cloud tracked in each 
block of clouds. 



Simulation 



m2chi2MR4 


0.55 


± 


0.07 


m3chi2MR0.25 


0.92 


± 


0.07 


m3chi2MRl 


0.83 


± 


0.14 


m3chi2MR4 


0.62 


± 


0.07 


m3chi2MR4_R16 


0.72 


± 


0.02 


ml0chi2MRl 


0.86 


± 


0.22 


m!0chi2MR4 


0.62 


± 


0.06 



that the destruction of the clouds leads to the highest den- 
sities in the remnant occuring over the outer half radius 
(in contrast, when there is no mass loading, a thin dense 
shell forms at the forward shock). The s e find ings have since 
been supported bv lDvson fc Hartquistl || 19871 ). who reported 
a similar "thick shell" morphology in their similarity solu- 
tio ns, and by th e addit ional numerical simulatio ns presented 
bv lDvson et ail < |2002f l and lPittard et all (|200ot ). The X-ray 
emission in these cases becomes softer and more extended. In 
other work, I Arthur fc Hennevl (|l996l ) studied the effects of 
mass loading by hydrodynamic ablation on supernova rem- 
nants evolving inside cavities evacuated by the stellar winds 
of the progenitor stars. They showed that the extra mass 
injected by embedded clumps was capable of producing the 
excess soft X-ray emission seen in some bubbles in the Large 
Magellanic Cloud. We conclude, therefore, that cloud de- 
struction by ablation, as in the simulations presented in our 



paper, can be looked for by searching for its affect on the 
X-ray emission of supernova remnants. 

The interaction between a wind and individual clouds 
is also seen in observations of planetary nebula. For in- 
stance, in NGC 7293 (the Helix nebula) long molecular tails 
and bright cresent-rimmed clouds are spectacularly resolve d 
l |0 'Dell et al.ll2005l ; iHora et alfood : iMatsuura etal]|2007l ) . 
The tails in the outer part of the nebula are less clear due 
to projection effects, but appear to be a separate population 
displaying wider opening angles. This may reflect changes in 
the diffuse flow past the clouds, possibly due to the mate- 
rial stripped off the clouds further upstream. Observations 
indicate that the flow past the clouds is mildly supersonic. 
Although numerical simulati ons are able to ma tch the ba- 
sic morphology of the tails (|Dvson et al.ll2006T ). dedicated 
3D simulations of clouds with very high density contrasts to 
the ambient flow are required in order for further insight to 
be gained. 

On larger scales, we note that there is substan- 
tial support for mass-loading in starburst regions. Broad 
emission-line wings are seen in many young star form- 
ing regions in c ludin g 30 Doradus ( Chu fc K ermicutt 1 19941 ; 
iMelnick et"afl I1999TI NG C 604 (lYang et al.l Il996l ). and 



19941) . 



NGC 2363 |Rov et aljfl992l ; iGonzalez-Delgado et al, 
and more distant dwarf galaxies (e.g. iMarlowe et al 
Izotov et al. 1 1 19961 ; iHomeier fc Gallagher! 1 19991; I Sidoli et al.l 



1995; 



20061 ). More recently. IWestmoquette et all (|2007al lbl) re 



ported on the broad emission-line component in the dwarf 
irregular starburst galaxy NGC 1569. Although the nature 
of the broad lines has yet to be fully determined, evi- 
dence is mounting that it is associated with the impact 
of cluster winds on cool gas knots. It therefore traces 
both mass-loading of wind material and mass entrainment. 
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Further support for mass-loading comes from the analy- 
sis of hard X-ray line emiss ion within starburst regions. 
IStrickland fc Heckmanl (|2009h determine that as much ma- 
terial is mass-loaded into the central starburst region of M82 
as is expelled by the winds and supernovae which pressur- 
ize the region. A key future goal is the development of nu- 
merical simulations of the multi-phase gas within starburst 
regions, and predictions for broad emission-line wings and 
X-ray emission. 

Starbursts are often associated with galactic out- 
flows. These flows are observed to be filamentary. The 
consensus view is that the clumps at the heads of the 
filaments are material which is ripped out of the galac- 
tic disk as the galactic wind develops (i.e. represent- 
ing additional, distributed, mass- loading) . In some cases 
the filaments appear to be confined to the edges of 
the outflow ijShopbell fc Bland-Hawthornl Il998h . while in 
other objects they appear to fill the interior of the wind 
[jVeilleux fc Bland-Hawthornl fl997l ) . It is clear that mate- 
rial stripped from the Ha emitting clouds is entrained into 
the outflow (see, e.g. ICooper et al.l [20081 ). but the exact 
amou nt is notoriously difficult to measure ()Veilleux et al] 
120051 ). Again, future dedicated simulations are needed to 
help address this issue. 

A key question concerns the ultimate fate of gas 
within galactic outflows. Observations and simulations in- 
dicate that the majority of the energy in galactic winds is 
in the kinetic energy of the hot gas, while the mass in the 
outflow is dominated by the warm photoionized gas. The 
former has a good chance of escaping the gravitational po- 
tential of the host galaxy, while the latter in many cases 
is unlikely to do so. Even when simple energetic arguments 
suggest that the ISM can be completely expelled from a 
starbursting galaxy, whether it will actually occur depends 
strongly on the geometry and m ultiphase nature of the ISM 
(see, e.g.. iHeckman et al.lll995l ). For instance, if a centrally 
concentrated starburst occurs in a galaxy with a disk-like 
ISM, blowout of the superbubble along the minor axis can 
allow the bulk of the ISM in the disk to be retained by 
the galaxy. With a multiphase ISM the diffuse intercloud 
medium may be ejected while large dense clouds remain in 
the disk. 



Recent work by lOppenheimer fc Pavel (|2008l ) indi- 
cates that the range of galactic winds is primarily deter- 
mined by the interaction of the wind with the ambient envi- 
ronment, with the gravity of the parent galaxy playing a less 
significant role. Their simulations also show that across cos- 
mic time the average wind particle has participated in a wind 
several times. However, further investigations are needed, 
since their simulations currently lack the resolution required 
to make accurate quantitative predictions of the slowing of 
the winds and wind recycling. Furthermore, the wind's mul- 
tiphase nature must be addressed. Studies like ours are rel- 
evant to this work since the destruction and acceleration 
timescales which we find for our clumpy region have some 
bearing on the mixing and stalling timescales of a galactic 
wind. Such simulations could be test ed against observations 
of the extent of g alactic winds (e.g. iTumlinson et al.l [20T1I : 
iTripp et alj|201lK . 



5 CONCLUSIONS 

We have performed a detailed investigation of a shock run- 
ning through a clumpy region. We find the following key 
behaviour: 

• The stripping of material from the clouds "mass-loads" 
the post-shock flow and leads to the formation of a dense 
shell. Fully-mixed material within the shell reaches a maxi- 
mum density, after which the shell grows in width. The shell 
expands and its density drops once there are no more clouds 
to mix in. 

• The evolution of the shock can be split into several dis- 
tinct stages. During the first stage, the shock decelerates. 
Then in some cases its velocity becomes nearly constant. Af- 
ter deceleration and the constant speed phase, if it occurs, 
the shock accelerates and finally approaches the speed it had 
in the uniform medium before it encountered the clumpy 
region. A steady stage does not always occur (e.g., if the 
clumpy region is not very deep). 

• When the mass-loading is sufficient, the flow can be 
slowed to the point that the shock degenerates into a wave. 

• The clumpy region becomes more porous as the number 
density of clouds is reduced, which occurs for lower values 
of MR, and/or increased values of X- 

• A great deal of turbulence is generated in the post- 
shock flow as it sweeps through the clumpy region. Clouds 
exposed to this turbulence can be destroyed in only 60% of 
the time needed to destroy a similar cloud in an "isolated" 
environment. The lifetime of downstream clouds decreases 
with increasing MR. 

We have determined the necessary conditions (in 
terms of the cloud density contrast and the ratio of cloud 
mass to intercloud mass) for a clumpy region to have a sig- 
nificant effect on a dif fuse flow. The l ifetim e of clouds is a key 
factor in this respect. IPittard et al] (|2009l ) first showed that 
clouds can be destroyed more quickly when overrun by a 
highly turbulent environment, although the strength of this 
turbulence was treated as a free parameter. In contrast, in 
this paper we have presented the first self-consistent simu- 
lations of a highly turbulent flow overrunning clouds, where 
the enhanced turbulence is a natural consequence of the flow 
overruning clouds further upstream. Since our simulations 
are purely hydrodynamic, the next step will be to investi- 
gate the effects of magnetic fields and thermal conduction on 
the flow dynamics. In future papers we will also investigate 
the ability of a flow to force its way through a finite-sized 
clumpy medium, and determine how this depends on the 
ratio of the mass injection rate from the clouds to the mass 
flux in the wind. 
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